Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 May 2012 (MN lAT^ style file v2.2) 



Introducing PHAEDRA: a new spectral code for 
simulations of relativistic magnetospheres 

Kyle Parfrey^'^*, Andrei M. Beloborodov^, and Lam Hui^'^ 

^Astronomy Department, Columbia University, 550 West 120th Street, New York, NY 10027, USA 

^ Institute for Strings, Cosmology, and Astroparticle Physics (ISCAP), Columbia University, New York, NY 10027, USA 

^Physics Department and Columbia Astrophysics Laboratory, Columbia University, 538 West 120th Street, New York, NY 10027, USA 



2 May 2012 



ABSTRACT 

We describe a new scheme for evolving the equations of force-free electrodynamics, 
the vanishing-inertia limit of magnetohydrodynamics. This pseudospectral code uses 
global orthogonal basis function expansions to take accurate spatial derivatives, allow- 
ing the use of an unstaggered mesh and the complete force-free current density. The 
method has low numerical dissipation and diffusion outside of singular current sheets. 
We present a range of one- and two-dimensional tests, and demonstrate convergence 
to both smooth and discontinuous analytic solutions. As a first application, we revisit 
the aligned rotator problem, obtaining a steady solution with resistivity localised in 
the equatorial current sheet outside the light cylinder. 

Key words: magnetic fields — plasmas — methods: numerical — relativistic processes — 
pulsars: general — MHD 



1 INTRODUCTION 

We present a code for simulations of force-free electrody- 
namics: PHAEDRA (Pseudospectral High-Accuracy Electro- 
Dynamics for Relativistic Astrophysics). The systems we 
study are 'force-free' in the sense that the Lorentz force 
density vanishes everywhere, because the electromagnetic 
fields are strong enough that hydrodynamic forces can be 
neglected, resulting in self-balancing electromagnetic fields. 
Relativistic force-free electrodynamics has long been recog- 
nised as the appropriate limit for describing the magneto- 
spheres of neutron stars and black holes, yet only recently 
has a concerted effort begun to study it with direct numer- 
ical simulation. 

This infinite-magnetisation, or vanishing-inertia, limit is 



the appropriate one for the magnetospheres of pulsars ( Gol-| 
dreich fc Julian|1969||Contopoulos, Kazanas, fc Fendt|1999[ 



Gruzinov 2006 Spitkovsky 20061, and magnetars, whose 



persistent and transient high-energy emission may be due 
to the distortion, reconnection, and dissipation of force- 
free fields ([Thompson fc Duncan||l995[ |Lyutikov|[2006l |Be^ 
|loborodov|2009[ ). Force- free electrodynamics is the standard 
tool to study the extraction of rotational energy from black 



holes ( [Blandford fc Znajek] |1977| [MacDonald fc Thorne 
|1982[ ), where the magnetic field is thought to be supplied 
by a conducting accretion disc. The natural self-coUimation 
of electromagnetic fields make them attractive candidates 
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for explaining relativistic jets in quasars and active galactic 
nuclei, whose high Lorentz factors suggest low baryon load- 
ing and electromagnetic dominance (Blandford 1976); simi- 
lar Poynting jets may be responsible for gamma-ray bursts 
(Meszaros & Rees 19971. An argument can be made that 



all ultra-relativistic outfiows are essentially electromagnetic, 
rather than gas dynamical ( Blandfo"rdl|2002 1 . 



In any case, it is clear that there exists a wealth of 
challenging problems for the field of electrodynamic numer- 
ical simulation. Direct time-dependent simulation is valu- 
able because it permits the study of general realistic initial- 
value problems, without the restrictions, like self-similarity 
or stationarity, that are often necessary in analytical models, 
and because it naturally tests the stability of field configura- 
tions, a question often unanswered by steady-state numeri- 
cal work. 

Several time-dependent force-free electrodynamics 
codes exist, both those using finite differences (Spitkovsky 
[2006; Kal apotharakos fc Contopoulos|2009||Palenzuela et al 



approach (|Komissarov||2004a 


|Cho||2005 


|Asano, Uchida, &| 


Matsumoto||2005 McKinney 


2006a Yu 


20111. Our numer- 



ical scheme is entirely different, and complementary, being 
based on orthogonal basis function expansions. 

Previous codes have large numerical dissipation or dif- 
fusion, introduced either because they do not maintain 
E ■ B = self-consistently or through the intrinsic diffusiv- 
ity of the method, while force-free problems often demand 
long simulations, as the fields may evolve over many wave- 
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crossing times. It is desirable to have a method which can 
be run for long times without intrinsic dissipation, captures 
discontinuities, and accurately describes fast dynamics. 

The crucial question one asks of a force-free configura- 
tion is that of its stability, the onset of instability commonly 
leading to a dramatic rearrangement of a magnetosphere, 
sometimes involving explosive reconnection. Spectral calcu- 
lations tend to have less numerical noise than those of com- 
parable finite-difference or finite-volume ('local') schemes; 
this noise can erroneously trigger instability. In a study of 
Sweet-Parker reconnection, the spectral magnetohydrody- 
namics (MHD) code is found to be largely immune to the 
secondary island formation, caused by a tearing-mode in- 
stability, that is found using local methods for the same 
problem ( |Ng fc Ragunathan|20TT| ). 

In this paper, we describe a code for axisymmetric sim- 
ulations, in flat space-time. It has been designed in such a 
way as to be extensible with minimal restructuring to a fully 
three-dimensional setting, in curved space-time. 

Force-free electrodynamics, and its relation to relativis- 
tic MHD, is discussed in Section [2] Section |3] contains a de- 
tailed description of the code and its practical implementa- 
tion, including some background on each of its components. 
A range of one- and two-dimensional test problems are pre- 
sented in Section |4] including convergence tests in realistic 
scenarios. The aligned rotator is examined in Section [5] Fi- 
nally, we discuss our results in Section [6j and outline some 
promising areas of future research. 

Note that in Section |3] we distinguish between the con- 
travariant, , and covariant, Fi, components of a vector 
field, while when we discuss results in Sections [4] and [5] we 
refer only to the components in an orthonormal basis, also 
written Fi. 



2 FORCE-FREE ELECTRODYNAMICS 

The system of force-free electrodynamics is the vanishing- 
inertia, or, equivalently, ultra-relativistic, limit of magneto- 
hydrodynamics. The latter can be written as 



Vp*^""' =0, 



VpJ — ^ y-l- (m) + -'(f) 



0, 



(1) 



(2) 



where *F'^^ is the Maxwell tensor (Hodge dual of the 
Faraday tensor F^"), and T^^j and T^^^'^ are the energy- 
momentum tensors of the matter and electromagnetic fields, 
respectively. If the matter contribution to T'^" can be ne- 
glected, equation ([2| simplifies to 

(3) 



(f) ^ 



Combining 



Tr.: = F'^^F-"^^^[F^^F"^)9^-', 



(4) 



where g^j^ is the metric tensor of space-time, with the inho- 
mogeneous Maxwell equations. 



V.F"" = J", 
one finds that equation (|3| becomes 

F^^r = 0, 



(5) 
(6) 



which states that the Lorentz force density vanishes (e.g. 
Komissarov'2002) . Equation ([6| can alternatively be derived 



by postulating, in a frame in which the electric and magnetic 
fields are parallel, that the electric field vanishes and the 
current flows along the magnetic fleld. 

Standard MHD codes can become inaccurate, and even 
crash, when the plasma's magnetisation is large, because the 
numerical error in the electromagnetic energy density be- 
comes comparable to the matter energy density ( [Gammie^ 
McKinney, fc T6th||2003[ |Komissarov|2004b| |. They also re- 
quire the (usually poorly constrained) matter distribution 
to be set at the beginning of the simulation, and maintained 
throughout, sometimes by ad-hoc matter injection. Force- 
free codes do not experience these difficulties. 

We move now to a 3-1-1 space-time point of view. Equa- 
tion ([gJi becomes 



PcE + J X B = 0, 



(7) 



where E and B are the electric and magnetic fields, and J 
the current density. We can see that 



E- B = 0; 



(8) 



the electric field is 'degenerate.' This condition, together 
with V ■ B = 0, implies that the system of electromagnetic 
fields has only four independent components. Likewise, 



E-J = 0- 



(9) 



there is no Joule heating, the system is dissipationless, and 
formally (conditionally) hyperbolic. 

The evolutionary Maxwell equations, equations ^ and 

are the familiar 



dtB = -V X B, 
dtE = \/ X B - J, 



(10) 



using Heaviside-Lorentz units with c — 1; this just means 
that our current density is 4n times the current density 
in Gaussian units. In MHD, an additional relation must 
be given for the current, closing the equations. In force- 
free electrodynamics, the current is uniquely determined 
by equations ^ and ( 10 1 , together with the condition 
dt {E-B) = 0, to be ( |Gruzinov|1999 l 



j^ B.VxB^E-VxE ^^^^ExB ^^^^ 
Ohm's law in force-free electrodynamics is therefore essen- 



tially geometrical. The first term in equation (111, is the 



conduction current parallel to B, which maintains the de- 
generacy condition, equation ([sjl. The second term is the 
drift current, being in the form pei'drift, where 

Ex B 



Vdrift 



B2 



(12) 



is the velocity of the magnetic field lines. It is apparent that 
there is a second condition. 



B-E^> 0, 



(13) 



equivalent to requiring the drift velocity to be less than the 
speed of light; since charged particles cannot cross field lines, 
this is a requirement if we assume that a macroscopic matter 
velocity field exists. Equations ([sjl and ( 13 1 are commonly 
referred to as the 'force-free conditions.' The second con- 
dition implies that the electromagnetic fleld is intrinsically 



© 0000 RAS, MNRAS 000, 000-000 



PHAEDRA; a new spectral code 3 



magnetic, in that a frame exists in which the electric field 
vanishes. 

It is possible for the fields to evolve from a configu- 
ration in which this second force-free condition is satisfied 
everywhere to one in which it is violated at some point, line, 
surface, or volume. This local breakdown of the force-free 
approximation is necessarily accompanied by dissipation, as 
degeneracy is broken and E ■ B ^ 0. In general, any con- 
figuration having field lines of different topology, such as 
the open and closed lines of the pulsar magnetosphere, will 
have points, lines, or surfaces at which \B\ = 0, violating 
the second condition ( |Uchida|1997| ). These sites of force-free 
breakdown are especially interesting, as the localised dissi- 
pation may be responsible for observed radiation. 

Force-free electrodynamics supports two classes of 
waves, which we describe in a frame in which E = (Pun 



sly|2003[ ). Fast waves, equivalent to vacuum electromagnetic 
waves, propagate isotropically with both phase and group 
velocities equal to the speed of light. They do not carry 
any charges or currents. Alfven waves can carry charges and 
currents, have phase velocity «phase = ±ccos^, where 9 is 
the angle between B and the wave vector, and have group 
velocity equal to c and directed along B, iJgroup = ±cB /B. 

The 'field-evolution' approach we take is not the only 
way to write the evolutionary equations of force-free electro- 
dynamics. One could equivalently evolve the drift velocity 
or the Poynting flux vector, S — E x B, instead of the elec- 
tric field. The equations can also be written, using Euler 
potentials, as a Hamiltonian system ( Uchida|1997 l, or in an 
axionic formulation ( [Thompson fc Blaes|1998| ^ 



3 NUMERICAL METHOD 

3.1 The pseudospectral method 

A function u{x,t), the solution of a time-dependent partial 
differential equation, can be expanded in terms of a set of 
orthogonal spatial basis functions (f)k[x), 



UN 



(14) 



where ak(t) are time-dependent expansion coefficients; ujv 
is an approximation to the function u for some choice of 
basis functions and truncation TV. Spatial derivatives can 
be taken by analytically differentiating -ujv, since the exact 
derivatives of the basis functions are known. Considering 
first an equation in x only, 

Du{x) = f(x), 

where D is a general differential operator and / is a forcing 
function, we can think of solving this equation by minimis- 
ing the residual R: R = Dun — f ■ In what sense we choose to 
minimise R will determine the kind of spectral method we 
construct. In the Galerkin (sometimes just called 'the spec- 
tral') method, the residual is made orthogonal to the basis 
functions: {4>k, R) = 0, k = 0, . . . , N — 1, where the brackets 
indicate an inner product, {f,g) = / u}{x)f{x)g(x)dx, over a 
weight function u}{x). Since the first A'^ spectral coefficients 
are exact, Wiv can be considered to be a truncation of the 
infinite series expansion. 

In the pseudospectral method, the residual is made zero 



at a set of 'collocation' points, {xi}: {S[x ~ Xi], R) — 0, i — 
0, . . . , N — 1, where 5{x) is the Dirac delta-function. The 
resulting is then an mterpolant of the true it, at the 
chosen grid points. It can be shown that, if the collocation 
points are chosen as the abscissas of a Gaussian quadrature 
associated with the basis set and this quadrature rule is 
used to calculate the inner products, then the Galerkin and 
pseudospectral methods are equivalent for linear problems; 
the error penalty for choosing interpolation over truncation 
is at worst a factor of two, for trigonometric functions and 
Chebyshev polynomials ( Boyd|[2"001 1. 

Gaussian quadrature of a function / over a weight func- 
tion UJ, 



f{x)uj{x)dx = ^ Wif{xi) 



(15) 



is accomplished by finding the corresponding set of N 
weights {wi} and N interpolation points {xi}; the pay off 
for being restricted to these interpolation points is that the 
resulting formula is exact for all f{x) which are polynomials 
of degree 2A — 1 or less. The weight function determines 
the basis functions; for example, the Chebyshev polynomi- 
als are those polynomials which are orthogonal with respect 
to the weight uj{x) — 1/Vl — x'^ on the interval [—1, 1]. For 
periodic f{x), the composite-trapezoidal and midpoint rules 
are Gaussian quadratures with an equispaced grid, the cor- 
responding basis being trigonometric functions. 

The pseudospectral method can also be thought of as 
the limiting case of increasing-order finite-difference meth- 
ods, where the derivative stencil now extends over all grid 
points. In particular, the Fourier pseudospectral method on 
a periodic uniform grid is recovered by a finite-difference for- 
malism as stencil width (and formal order of accuracy) goes 
to infinity ( Fornberg|1996 \ . This approach gives a dense dif- 
ferentiation matrix, whose application requires 0{N^) op- 
erations. However, identical derivatives can be calculated 
for the Fourier and Chebyshev basis sets by using the fast 
Fourier transform (FFT), which requires only 0{N In N) 
operations; this is sometimes referred to as the 'trans- 



form method' (Orszag 19701. In this forward FFT 



gives the expansion coefficients, {ak}, from which can be 
constructed the coefficients for the derivative series {a'k}: 
u'n{^) ~ X^fcLo '^fc'^'=(-^)- Fourier basis, (j)k{x) = e'*"^ 

and so = ikak ; with Chebyshev polynomials a three-term 
recurrence relation is used. Finally the derivative at the grid 
points, u'pf, can be found with an inverse FFT. Given a func- 
tion to be differentiated, this procedure can be thought of 
as finding an interpolating function of order N, at a set of 
N points, and taking the exact analytic derivative of this 
interpolating function. 

The great benefit of these methods is that spectral 
approximation is exponentially convergent for sufficiently 
smooth functions: the error decreases faster than any power 
of the truncation N. It has generally been found that this 
carries over into exponential convergence of spectral solu- 
tions of PDEs, even those with fixed-order time marching. 
This accuracy has made spectral methods popular in many 
areas of physics, including meteorology, seismology, shock 
waves, and reactive flows. Astrophysical applications in- 



clude accretio n disc magnetohydrodynamics ( Chan, Psaltis, 
& Ozel 2005 2009| ) and general relativity (e.g. Gourgoul- 
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hon||1991| [Bonazzola et al.||1999| [Kidder et al.||2000| |Dim-| 



melmeier et al.|200 5 : Grand clement fc Novak|2009 1 . In engi- 
neering electrodynamics, the 'pseudospectral time-domain' 



method was introduced by Liu ( 19971, where it was shown to 



have much lower diffusion and dispersion error than finite- 
difference methods, and to require either two (Fourier) or tt 
(Chebyshev) points per wavelength for adequate resolution, 
in comparison to eight to sixteen points for finite differences. 
These lower required grid densities make a spectral method 
more efficient for achieving a given accuracy, despite the 
higher number of operations per grid point. 

3.2 Spatial discretisation 

In order to simplify eventual extension to curved space- 
time, we adopt a scheme which allows the use of an arbi- 
trary spatial metric. We store, and advance in time, the 
contravariant vector components of B and E\ curls are 
taken with {VxPy = {1/ ^)e'^''djFu and diver gences with 
V ■ F = (l/y7)9i(^/7-f ), where 7 is the spatial metric de- 
terminant, e*-'* the Levi-Civita symbol, and F stands for B 
or E. The quantity to be differentiated, either Fk or , 
is first calculated at each point from the contravariant com- 
ponents F^ and the metric, then expanded in orthogonal 
basis functions. This method requires more forward trans- 
forms than one where the derivatives are simplified by the 
chain rule, since for example both Fr and ^F^ must be 
transformed into spectral space, but we find it to be more 
accurate. 

Our grid is defined in axisymmetric spherical coordi- 
nates {r,6), with A'^ collocation points in r and L points in 
9. We will use i and j to index grid points, and n and I to in- 
dex wavenumbers, along each direction; i, n = 0, . . . , A'^ — 1, 
and j, Z = 0, . . . , L — 1. 

In the radial direction we use Chebyshev polynomi- 
als, r„, and the collocation points are chosen to be the 
Chebyshev-Gauss-Lobatto nodes, 



i = Q,...,N-l, xG[-l,l]. 

(16) 

These can be mapped directly onto the physical coordinate, 
ri = Tin + (rout - nn)(l + Xi)/2, r G [nn, rout], or via an ad- 



N - 



ditional coordinate mapping (section 3.31. Since the Cheby- 
shev polynomials are mapped cosine functions. 



T„(cos[g]) = cos(nq) , 



(17) 



with this choice of grid the Chebyshev transform of a func- 
tion / can be performed with a fast cosine transform: 



N-1 

= ^ /„ cos 



TV- 1 



(18) 



In the meridional direction we expand in sine or co- 
sine functions, depending on whether the vector component 
in question is an even or odd function across the pole (its 
parity). This is related to the 'double Fourier' method of ex- 
panding functions on a sphere ( Merilees|1974 Qrszag|1974 1 , 
which is attractive because it avoids the slow Legendre trans- 
form required by spherical harmonics. The following are 



even, and can be expanded in cosines: Fr, F^ , F^, F'*, 
whereas odd functions that can be expanded in sines are 
Fg, F^ , y/^F"^, y/^F"^. Here we are only concerned with ax- 
isymmetric modes; in general the parity will depend on 
whether the azimuthal wavenumber is even or odd. To avoid 
solving the equations directly on the poles we use a shifted 
grid: 



j + 1/2 . 
L 



J =0, 



1 



(19) 



For instance, the covariant radial component of the 
magnetic field, once formed by direct index lowering with 
the metric, can be expanded as 

JV-l L-l 

Br=Y,Y. '^r.iTr^ir) cos{W) , (20) 

n=0 i = 

and the combination y^E^ , required to calculate V ■ E, as 

N-l L-l 

^/7S'■ =Y.Y1 «"'^"(^) sin(Ze) . (21) 

n = 1 = 

Once the coefficients a„i of a function have been found, 
by taking the forward transforms in both directions, the co- 
efficients of the derivative series a'„i can be calculated. For 
differentiation by 9 this is simple: a'„i = —la„i for functions 
expanded in cosines, and a'„i — lani for those expanded in 
sines. Radial differentiation requires the three-term recur- 
rence relation, relating the Chebyshev coefficients of a func- 
tion to those of its derivative, 

= 

ajv-2,i = 2(Ai' — l)ajv-i,! 
c„_ial,_i^, = a^i+i,! + 2na„,; , (22) 
where Co = 2, and all other Cn, = 1. 

3.3 Coordinate mappings 

The standard Chebyshev nodes, Xn, are not a suitable radial 
grid for our problems. They are very strongly clustered near 
the endpoints, leading to a time step restriction for hyper- 
bolic problems that goes like At ~ 0{1/N'^), which makes 
obtaining high spatial resolution in the rest of the domain 
unnecessarily expensive. A coordinate map can be used to 
relieve the endpoint clustering, and also put more nodes in 
the part of the domain where higher resolution is required. 
Consider a map from x to a new coordinate y: 

y = gix) 

d_ _ 

dy g' dx ' 

where g' = dg/dx. The derivative coefficients can be calcu- 
lated as before, and the mapping of the derivative achieved 



by a multiplication in real space. The arcsine map (Kosloff 
fc Tal-Ezer| 19931 ), 



: S'asin(a;) 



arcsin(aa;) 
arcsin(Q:) 
a 



arcsm(Q) Vl - a^x^ 

where q is a constant between zero and one and x^sin G 
[—1, 1], will create an almost equispaced grid with less clus- 
tering at the endpoints; the grid stretching becomes more 
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pronounced as a is increased towards unity. Care must be 
taken in choosing a: since the map is singular at the end- 
points, choosing a too large would impair the convergence of 
the series. The maximum discrepancy between a functioij^ 
f{z) and its truncated series representation Pn{z) is 



:\f{z)-PN{z)\=ce. 



(24) 



where c is a constant which depends on / but is independent 
of A^, and e is related to a by 



1 - Vl - a'^ 



a = sech 



lne| 



(25) 



Therefore choosing e to be sufficiently small will usually 



make the singularity harmless (Don & Solomonoif 1995 
Mead fc Renaut|2003 1 . Equation ( 24 1 determines the largest 
pointwise error, which usually appears at, or close to, a 
boundary; the error can be smaller than e in the interior of 
the domain, which we confirm in our tests in Section |4] For 
fixed e the minimum grid spacing only decreases as 0(1/N). 

The nearly equispaced grid resulting from the arcsine 
map is still not perfect for computations in spherical coor- 
dinates when rout ^ rin, because the lines of constant 6 
converge towards the centre of the domain, and so too much 
radial resolution is used where the angular resolution is low, 
and not enough where it is high. We would prefer a grid 
where Ar ^ r AS over most of the domain. To construct 
this we use a combination of the arcsine map and a smooth 
algebraic stretching: 




9aig(a;) 



(26) 



Figure 1. Grid spacing versus radius: Chcbyshev nodes with a 
combination arcsine -|- algebraic mapping; N = 256, L = 155, 
Q = 0.6, € = 10-". 



an explicit Runge-Kutta integrator. We mainly use a third- 
order, three stage, method ( Fornber"g|[r996| ) . If the ODE to 
be solved is 



where Q is a constant and Xaig G [—1,1]. The grid is then 
linearly transformed to the desired physical coordinates 



du 
dt 



= fiu,t), 



-Rin)(l + a;aig), r G [Rin, Rout] ■ (27) 71 + 1, it"+\ from that at time n by 



r = Rin + - {Rou 

We generally use Q ~ 0.1 — 1, and set the map-induced error 
to e~ 10-^-10-^^ It appears to be preferable to have the 
radial grid spacing somewhat smaller than the meridional 
spacing throughout most of the domain. Figure [T] shows the 
inter- nodal spacings of an example grid. 

Some of our models also include a coordinate transfor- 
mation in the meridional direction, in order to increase the 
resolution around a reconnecting current sheet in the equa- 
torial plane. To this end we employ the 'Kepler-Burgers' 
mapping ( [B^p92l ), 



G {B^,E^}, then we can find the solution at time step 





= Atf{u",t") 






= Ai/(M" + ^fc(^',t" + 










n + 1 


= u" + -(fc<''+3fc'^>). 
4 





■ = e+|sin(2e). 



(28) 



where 9 is evenly spaced on [0, tt] and d is the new stretched 
coordinate. The constant 7 controls the degree of stretching; 
when using this map we set 7 ~ 0.3 — 0.5. 



3.4 Time evolution 

When the spatial derivatives have been found the system of 
equations becomes a set of ordinary differential equations in 
time, one for each vector component, which we solve with 



^ The complex function f{z) is the analytic continuation of the 
real function f{x). 



(29) 



Since fc'^-* is not required once fc*-"^' has been calculated, 
the latter can overwrite the former in memory, and so 
this method only requires storage for two intermediate ar- 
rays. Explicit time advancement is subject to the Courant- 
Friendrich-Lewey stability constraint on the time step, At 
CcFL'Jmin, whcrc Smin is the Smallest grid spacing, in our 
case always found at the inner surface. The factor Ccfl is, 
in general, problem and time-integrator dependent; we find 
that equation (291 is stable for Ccfl 5; 1, and hence always 

set At — 5min. 

We compared the above with a fifth-order, six stage, 
Runge-Kutta integrator. The improvement in accuracy is 
negligible, even for our stringent test with the Michel 
monopole solution, indicating that time-stepping errors are 
subdominant. This is unsurprising given that a fine grid in 
the radial direction is required near the stellar surface to pro- 
duce highly-accurate solutions, and so the stability-limited 
time step produces small errors. 
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3.5 Spectral filtering 

The use of spectral methods for non-linear hyperbolic prob- 
lems comes with two principal difficulties. The first is the 
build-up of power at high wavenumbers due to non-linear 
couplings between lower wavenumbers. Hyperbolic prob- 
lems have no explicit dissipation in the equations of motion, 
and spectral methods have very low intrinsic dissipation, so 
these high modes do not decay, and can lead to the break- 
down of the scheme by the aliasing instability. Aliasing oc- 
curs because a discrete transform of point values mistakes 
high-frequency components for low; all modes g"('^+9^)^^ for 
g = 0, ±7r, ±27r . . ., will be identical when represented on a 
grid of N points in the interval [— 7r,7r]. Non-linear terms 
produce high-frequency modes, which are then mistaken as 
low-frequency modes by the discrete transforms; these phan- 
tom low-frequency components are then similarly combined 
to generate more power at high frequencies, a clearly unsta- 
ble cycle. 

The second difficulty is the ability of non-linear cou- 
plings to create discontinuities in a solution which was pre- 
viously smooth. Spurious oscillations appear in the spectral 
interpolant when a function is not resolved by the gric|^ this 
is the Gibbs phenomenon. A jump discontinuity causes 0(1) 
errors in its immediate vicinity, and reduces the convergence 
rate to first order elsewhere. 

Both of these difficulties can be largely overcome with 
spectral filters. If a function u is expanded in a basis <j)n, its 
filtered approximation, J^u, is given by 



The use of the exponential filter to control aliasing was 



7V- 



(30) 



where m„ are the discrete expansion coefficients. IVandeven] 
(19911 showed that if the filter function, <j{ri), is unity at 
r] = 0, zero for all \ri\ ^ 1, and has at least 2p — 1 con- 
tinuous derivatives, then will converge with A'^ at 2p-th 
order even in the presence of a jump discontinuity, except 
very close to the jump. In addition, since it strongly damps 
the high modes, such a filter can prevent the onset of the 
aliasing instability if regularly applied to each field in a time- 
dependent simulation. We tried two variable-order examples, 
the erfc-log filter ( |Boyd|1996| , 

--II 1 



^eric-iog(r7) ^ ^erfc |2(2p)^/^r;^ ^^'^ | , (31) 



and the exponential filter ( |Majda et al.|1978 l 



0-exp(»7) 



(32) 



The latter can be made to fulfil approximately the require- 
ment of being zero for all I?)! 1 by setting a — om — 
— In(eM), where em is machine precision; we use ~ 35. 
We find the exponential filter to give more accurate results 
and to allow weaker filtering, and so use it exclusively. 



studied in detail by Hon & Li (20071, where they found 
that a high-order (2p = 36) filter can prevent instability in 
marginally-resolved fluid dynamics simulations while pro- 
ducing more accurate solutions than standard dealiasing 
methods. We also find this to be the case for our equa- 
tion system, even though the non-linear coupling is much 
stronger than in Euler's equations. For all science runs we 
use a filter with a = and 2p = 36, which appears to bal- 
ance well the confiicting demands of assuring stability while 
minimising unphysical dissipation. When very low numbers 
of modes are used (roughly A < 48 in any direction) the fil- 
tering order needs to be reduced somewhat. This high-order 
filter is applied to the coefficients of every derivative serie^ 
and directly to the coefficients of the solution itself at the 
end of each full Runge-Kutta time step. 

We turn now to the second issue, concerning stabil- 
ity and convergence in the presence of jump discontinuities. 
Current sheets, regions of formally infinite current density 
implying discontinuities in the magnetic field, are a generic 
feature of force-free electrodynamics. The danger is that the 
strong production of high-wavenumber power, or the non- 
linear interaction of the Gibbs oscillations with the solution, 
will lead to instability, or at least loss of convergence. Dissi- 
pation is required to prevent this, and ensure that the cor- 
rect entropy solution is selected; see [Gottlieb fc Hesthaven] 
( |2001[ ) for an extensive review of stability and convergence 
theory for non-linear hyperbolic problems. 

This dissipation can be effectively, and efficiently, pro- 



vided by a spectral filter (Tadmor 1993 Don 19941. For a 
time-dependent problem. 



du 



we can advance the solution one time step, and then apply 
a filter to the solution, u{t + At). Note that in this case 
we are using the term 'filter' more broadly than above; in 
particular, we will not require that ct(1) — 0. If we choose to 



use the exponential filter, equation (32 1, then the modified 



equation, taking into account the action of the filter, is 



if u is expanded in Fourier series, and 



du 

m 



= /(«) 



At A2p 



dx 



n 2p 



. + 0{At') , (34) 



if it is expanded in Chebyshev polynomials. The second re- 
sult follows from the relation 



dx 



T„{x) + n^T„(x) = . 



(35) 



In this case the dissipation decreases to zero at the bound- 
aries (recall x G [—1,1], and here u{x) is non-periodic), 
which is useful since then no additional boundary conditions 
are required ( |Boyd|1998| . 



Any set of points on a grid will be faithfully recovered follow- 
ing transforms into, and back out of, Fourier space; however, the 
interpolating function, constructed from the Fourier coefficients, 
can show spurious oscillations between the grid points, causing 
oscillations in the derivative at the grid points. 



When using Chebyshev polynomials it is important to filter 
the coefficients of the derivative series, a' ,, rather than those 
of the function before the recurrence relation is used, the a„; 
(Godon & Shaviv 1993|. For this reason our filtering operations 
are implemented inside the coefiicients-to-grid inverse transform. 
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Applying an exponential filter is therefore similar to 
adding a hyperviscous term to the equation, where the mag- 
nitude of the hyperviscosity is 

= aS^ ■ (3*^) 

The action of the filter is equivalent to an implicit time 
integration of the hyperviscous term, and so no additional 
stiffness is added to the equations. 

The similarity of the effect of the exponential filter to 
an explicit hyperviscous term allows us to import stabil- 
ity and convergence theory derived using such terms. The 
spectral viscosity (SV) method (Tadmor 1989 19901 uses 



second-order viscous regularisation (2p = 2), and conver- 
gence is obtained by excluding an increasing fraction of 
low-wavenumber modes from the viscous term as A'^ is in- 
creased. Extensions to higher order, 2p ^ 4, followed for 



schemes based on Fourier (Tadmor 19931 and Chebyshev 
|[Malll998^ expansions, known as the super spectral viscos- 
ity (SSV) methods. These can be proven to converge to the 
correct solution of a scalar conservation law for 

ps; O(lniV) . 

Convergence can't be guaranteed for a system of equations, 
although it has been proven that if the scheme converges, it 



converges to the correct entropy solution (Carpenter et al 



20031. Despite the lack of solid theoretical results for non- 



linear systems, experience has shown that the method can 
be stable and convergent; multidimensional examples in- 



clude shock- vortex interaction ( |Don|1994 Sun et al.|2006 1, 
and problems involving both shocks and combustion (Don 



k. Gottlieb 1998 Gottlieb & Gottliebl 2005 1 , where spec 



tral methods were found to perform well in comparison with 
high-order shock-capturing schemes. 

Using equation (37 1 as a guide, we find from equa- 



tion ( 36 1 that the filter amplitude should scale as a = 
Qssv = CNAt. The time step scales as At oc 1/N (Sec- 
tion 3.3 1, and so qssv should be roughly constant. The fil- 



ter order, 2p, should only increase slowly with N, and so we 
set it to a constant as well. Numerical experiments confirm 
that fixed Qssv and p, determined by low-resolution simu- 
lations, lead to stable and convergent results as resolution 
is increased. We find best results are obtained with 2p — 8, 
Qssv ~ 0.01 — 0.1, corresponding to a weak hyperdiffusion 
which decreases in strength with resolution like A~^. The 
SSV filters are applied to every component of the electric 
and magnetic fields at the end of each full Runge-Kutta 
time step. 

To summarise the filtering procedure, we apply a very 
high order filter, with a = 35 and 2p ~ 36, to the inverse 
transform of every derivative series. At the end of each full 
time step, we apply both the previous filter and one with 
a ~ 0.01 — 0.1 and 2p = 8 to the field variables themselves. 

3.6 Post-processing 

The SSV method is exponentially convergent in any error 
norm, if p increases linearly with N; however stability will 
often not permit this, and in any case 0(1) errors will remain 
near any discontinuities. It has long been argued that high- 
order methods retain enough information to reconstruct a 



highly accurate (and sometimes spectrally convergent) solu- 
tion, even if oscillations are present and the pointwise errors 
are large. The Gibbs oscillations are not noise, and they can 
be safely removed, without destroying the accuracy of the 
underlying solution, with a post-processing step after the 
simulation has finished. There are many ways to do this: 



for example, real-space filtering or mollification ( Gottlieb & 
Tadmor|1985[ ), one-sided filters ( |Cai, Gottlieb, fc Shu|1992[), 
spatially-varying spectral filters ( |Boyd|1996||Tadmor fc Tan- 
ner|2005 1, and reprojection into a Gibbs-complementary set 



of basis functions ( Gottlieb et al. 1992 Gelb & Tanner 2006 1 



The reprojection method is particularly popular, and many 
examples exist of successful one-dimensional reconstructions 
of oscillation-free solutions (e.g. Shu & Wong 1995 Sarra 
|2003l[MalFLl|[2006l ). 

The problem with all of the above methods is that they 
require the locations of the discontinuities to be known ac- 
curately, which is particularly difficult in more than one di- 
mension. We were unable to find a sufficiently robust means 
of determining the number of unresolved features and their 
locations, given that the physics generates, and our scheme 
is capable of resolving, oscillations with wavelengths close 
to the grid scale. 

We show here an example of post-processing applied 
to a solution with a known discontinuity — the equatorial 
current sheet in the aligned rotator problem (Section [5|. 
Figure |2] shows (a) the original SSV solution before post- 
processing, (b) after applying the optimal spatially-varying 
spectral filter ( Tanner|2006 1 : 



o'opt(fc, N,x) = e 

ak^d{x) 
2N 



[ftJVd(a;)J 

E 

n=0 



(38) 



where d{x) is the distance to the nearest discontinuity and 
K and a are constants, and finally (c) using the digital total 
v ariation (D TV) spatial filter ( |Rudin et aL]|1992[ [Biirgel fc 
|Sonar|2002| . The DTV filtered solution u to a noisy variable 



u is found by minimising the fitted total variation energy 



A 



WdTV = ^ jVlt^l + {up - Ufj) 



(39) 



IVu^l 



where /3 ranges over all points in the dataset, 7 denotes 
each point's neighbours, and A is related to the expected 
noise level. The minimisation can be done by linearised Ja- 
cobi iteration. This method does not require the locations 
of the discontinuities, has a natural multi-dimensional form, 
requires only an estimate of the size of the oscillations to be 
removed, and has been applied to Chebyshev-based spec- 
tral methods with good results ( |Sarra|[2006[ ). However, we 
find it to perform poorly if resolved physical high-frequency 
oscillations are present, even with an adaptive local noise 
estimate. 

We do not use any post-processing technique for our 
results in Sections |4] and |5] with the exception of Fig. [Ts] 
where we apply a two-dimensional DTV filter. Our inten- 
tion is simply to highlight that any Gibbs oscillations which 
are present do not destroy the accuracy of the scheme. The 
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(a) No post-processing 
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(b) Optimal filter 




0.5 

(c) Digital total variation 

Figure 2. Toroidal magnetic field with discontinuity (current 
sheet), L = 255: (a) no post-processing, SSV only; (b) with opti- 
mal filter, K = 0.05, a = I; (c) with DTV reconstruction. 



pointwise errors may be large near a jump, and we do not 
claim spectral convergence of the field values on the grid, 
but the overall evolution of the system is still consistent 
with high-order accuracy. 

3.7 Force- free current 

The current density function required to close Maxwell's 
equations with force-free dynamics can be written as 



J \\ + J drift , 

B V X B 



E -V X E . 



B2 



J drift —V-E 



E X B 

(40) 



The effect of is to maintain the force- free condition 
E ■ B — 0. To date, schemes that rely on staggered grids 
do not explicitly include this term, because it demands that 
all components of the fields and their curls be evaluated, by 
interpolation if necessary, at the locations where the electric 
field components are defined. Instead, the effect of the par- 
allel current is mimicked by resetting the values of E on the 
grid such that E^ = 0. Since the inherent accuracy of spec- 
tral derivatives frees us from the need for staggered grids, 
we can evaluate the full current function without interpo- 
lation, and so include both Jy and Jdrift in the equations 
of motion. We also manually set E^^ to zero at the end of 
each full time step; this has essentially no discernible effect 
in most of the domain since J\\ keeps E ■ B very close to 
zero anyway, but gives slightly cleaner evolution close to the 
inner boundary at lower grid resolution. Specifically, E is 
projected, parallel to B, into a plane perpendicular to B: 



E - 



E-{E-B) 



B2 



(41) 



Retention of the parallel current results in much lower dif- 
fusion error; see the twisted dipole test in Section [4. 2. 2| 

Force-free evolution can lead to configurations in which 
the second force-free condition, — > 0, is violated. For 
example, in the equatorial current sheet of the axisymmetric 
rotating dipole (Section [5| all components of the magnetic 
field are zero, and any electric field results in the violation of 
the condition. An ideal force-free configuration requires E = 
in the current sheet. This can be simulated by immediately 
reducing the magnitude of the electric field if it is larger than 
the magnitude of the magnetic field, such that — — 
at the end of the operation. At each Runge-Kutta substep, if 
the condition is violated, we shrink the electric field vector, 
leaving its direction unchanged: 

'>2 

(42) 

This removal of electric field acts like a small, highly lo- 
calised, source of dissipation in current sheets. Physically, 
the removed electromagnetic energy would be converted to 
thermal energy and radiation; in these simulations it is sim- 
ply lost from the system. 



3.8 Boundary conditions 

A system of hyperbolic equations generally requires suitable 
conditions to be provided at the boundaries of the compu- 
tational domain. These boundary conditions will depend on 
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the spatial geometry and the physical problem under in- 
vestigation. We often wish to simulate a volume of force- 
free plasma surrounding a neutron star, and so describe the 
boundary treatment for this case. 

In these simulations, the inner boundary, r — r-m, cor- 
responds to the stellar surface, while the outer boundary 
at Tout should as closely as possible behave as a membrane 
which perfectly transmits outgoing waves without generat- 
ing any incoming waves. The behavioural boundary condi- 
tions at the poles are satisfied automatically by the choice 
of either sines or cosines as the basis functions in the merid- 
ional direction, as described in Section 13.21 



3.8.1 Inner boundary 

For our purposes the star is a perfect rigid conductor, all or 
part of which can be rotated. Since the star is a perfect con- 
ductor, the electric field will be zero in the rotating frame, 
E' = 0. The fields in the lab frame are given by 

B = B', 

E^-{nxr)xB', (43) 

where fl is the local angular velocity vector, and therefore 

- {nxr)x B. (44) 

Rotation about the z-axis corresponds to the application of 
an induced poloidal electric field. 



Equation ( 44 1 provides the relationship between E and 



B infinitesimally below the surface, but our first grid points 
are in the force-free plasma infinitesimally above. The nor- 
mal component of the magnetic field, , and the tangen- 
tial components of the electric field are continuous across 
the surface, and therefore are known. The required bound- 
ary values are B"" = B'"(6»), E" = -QB^'sme, and E"^ = 0; 
these are strongly enforced at every Runge-Kutta substep. 
The other components must be allowed to evolve undis- 
turbed, since they depend on unprescribed surface currents 
and charges on the star. 

A complication introduced by the combination of the 
above boundary conditions and the SSV filtering of the fields 
is the anomalous leakage of energy into the domain through 
the inner boundary. Consider the B^ field: at the end of 
each time step it is filtered, or smoothed, leading to a slight 
broadening of the -B' (r) profile, and hence diffusion of field 
from the boundary into the domain. According to the above 
boundary condition, stating that B^ is constant in time, this 
field at the boundary is immediately replenished, and so over 
time this cycle can increase the volume-integrated energy on 
the grid. 

A solution is to subtract the initial vacuum fields 
{Bo, Eq) from the dynamical fields (B, E), and evolve B — 
B — Bo and E = E — Eo. The current density and elec- 
tric field in the initial configuration are zero, and so, writing 
B = Bo + B, E = E, we get the equations 



dtB = -V xE 

dtE ^ V X B - J{Bo + B, E) 



(45) 



Note also that V x Bo = 0, and so VxB = VxB in 
the current density function — no additional derivatives are 
required. There is no requirement that these new variables 
be small, and in much of the magnetosphere they will be 



larger than the corresponding initial field. The initial field 
is never differentiated or filtered, there is no field leakage 
by the above mechanism since B^ = is the new magnetic 
boundary condition, and the energy conservation of the code 



is improved dramatically (see Section 4.2.21 



To reduce notational clutter, these variables will be nei- 
ther distinguished nor discussed elsewhere in this paper; the 
only addition they require is that the evolved magnetic field 
be temporarily added to a stored initial field before being 
passed to the current density function (which takes as argu- 
ments the electric and magnetic fields, and their curls). 

Finally, a small, slowly growing, anomalous toroidal 
magnetic field was found in certain simulations using a 
dipole base field, in the equatorial region immediately next 
to the star. This field appeared only when that region was 
'dead' (V x B = 0). We have attributed this feature to 
Alfven waves on under-resolved field lines, propagated with 
a numerical scheme with very low diffusivity. 

This phenomenon is completely eradicated by setting 
a small circular region of the poloidal plane, centred on 
{r,9) = (rin,7r/2) and with radius rdrift, to only use the 
drift current contribution to J in the equations of motion. 
Neglecting Jy makes the scheme in the 'drift' region suf- 
ficiently diffusive that the anomalous feature does not de- 
velop. Inside this region the electric field is projected to be 



perpendicular to the magnetic field, using equation (41 1, at 



the end of each full time step, exactly as in the rest of the 
domain. The radius of the drift region can be decreased with 
increasing resolution; we use rdrift ~ 0.25 rin. We have found 
that modifying the scheme in this small fraction of the do- 
main does not affect the solution in the rest of the domain. 

3.8.2 Outer boundary 

We implement the non-reflecting boundary condition at the 
outer boundary using an approximate characteristic decom- 
position. Spectral methods are very sensitive to 'incorrect' 
boundary conditions, and the simple zero-gradient condition 
that works well in low-order methods leads to instability. 



Characteristic-based boundary treatments (e.g. Abarbanel 
et al.|1991| |Godon|1996"| specify the out going characteristic 
variables using the calculated data on the grid, and set the 
incoming characteristic variables to zero. 

If we construct a six-component vector of the fields like 
q — {B,E}, then the one-dimensional equations of motion 
in Cartesian coordinates can be written in the general form 

dtq + Aa,q = 0, (46) 

where the matrix A (the flux Jacobian) should include con- 
tributions from the linear and non-linear terms. This matrix 
can be decomposed into its eigenvalues and eigenvectors, 

A = SAS~\ (47) 

A being a diagonal matrix of the eigenvalues, S containing 
the right eigenvectors, and the left eigenvectors. The 
characteristic variables, w, are found using the left eigen- 
vectors, 



c-l 
S q, 



(48) 



using which the equations of motion decouple, since A is 
diagonal: 



dtw + Adxtv = 0. 



(49) 
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In this one-dimensional case each component of w will move 
in either the positive or negative x direction, depending on 
the sign of the corresponding eigenvalue. The incoming vari- 
ables can be identified, set to zero, and the primitive vari- 
ables recovered using q — Sw. 

Rather than using the exact characteristic variables, we 
have implemented an approximate boundary condition using 
the characteristics of the vacuum Maxwell's equations. 

At each point on the outer boundary, construct a lo- 
cal Cartesian vector basis, with x in the radial direction, y 
along 6, and z along t^. We will then identify By = r_B*, 
Bz = rsmOB'^ etc., where the Cartesian components are in 
a normalised basis, and the spherical components are con- 
travariant as usual. The four propagating characteristic vari- 
ables are 



\wa/ 



\Ey 



B 

By 

B,J 



(50) 



The variables wi and W2 propagate in the negative r direc- 
tion, these are incoming variables and will be set to zero; 
U13 and propagate in the positive r direction and will be 
calculated from the field values on the grid. Inverting equa- 



tion ( 50 1, zeroing the incoming characteristics, and replacing 



Cartesian with spherical components gives 



(B^\ 

B"" 
E" 

\eV 



( ^"^3 \ 

U14 / sin Q 
\wa/ sin 9/ 



(51) 



the indicated fields at the boundary are replaced with these 
values at every Runge-Kutta substep, the radial components 
being left unchanged. 

This approximate boundary condition is stable and 
works well for predominantly radial waves, but can gener- 
ate artefacts when sizeable tangential waves are present. To 
prevent these from appearing we use a thin sponge layer 
next to the outer boundary, which absorbs outgoing waves 
( If et al.|1987l ). Any waves reaching the boundary are much 
attenuated, and so the vacuum boundary condition is a bet- 
ter approximation. Likewise, the approximate non-reflecting 
condition allows the use of a thinner and weaker absorbing 
layer than would be sufficient with a refiecting boundary. 

The sponge layer is introduced by adding a frictional 
term to Maxwell's equations, which become 



dtB = -V X f;-(j,(r)Bang 
dtE = V X B -J-as(r)E, 



(52) 



with the values cro ~ 0.5 — 1, 7 ~ 6, /3 ~ 4. This boundary 
treatment is robust, effective, and insensitive to the values 
of the sponge layer coefficients. 



3.9 Magnetic field divergence 

Maxwell's equations comprise two evolution equations, for 
the electric and magnetic field vectors, and two constraint 
equations. One of these constraints, Gauss's law, is automat- 
ically satisfied, since the charge density has been replaced by 
V ■ E in the drift current term in equation ( 40 1 . The other 



constraint is the solenoidal condition on the magnetic field: 

V - B = 0. 

In theory this should not be a concern, since the evolu- 
tion equation for B implies 9t(V ■ B) = —V ■ {V x E) = 0. 
However, in numerical schemes the operators for calculating 
divergences and curls usually do not satisfy V ■ (V x V) — 
exactly for any vector V, raising the worry that this 
truncation-level magnetic field divergence might build up 
over time. The presence of such magnetic charges can lead 
to unphysical forces and even instability. 

To maintain stability many MHD codes use constrained 
transport, which can ensure that some representation of 

V ■ i? is kept at machine zero (see T6th| |2000) for a re- 
view). These methods are incompatible with our global spa- 
tial derivatives. 

Another option is to evolve the magnetic field using a 
vector potential: writing B = V x A, the first evolution 
equation becomes dtA — —E; see Chan et al. ( 2009| for a 
spectral implementation. Although V • (V x A) is only zero 
to truncation error, this error would not grow over time. 
The disadvantage of a vector potential is the introduction 
of second-order spatial derivatives, and, in our case, an in- 
crease in the number of derivatives that must be taken. More 
problematically, we found this method to be less stable than 
the direct magnetic field-evolving method, especially at the 
boundary. 

It appears that the best approach in the context of a 
spectral method may be to do nothing, and rely on the accu- 
racy of the derivatives, and hence the smallness of the trun- 
cation error, to maintain the solenoidal condition to high 
precision. Let us define a normalised magnetic divergence, 



(V-B)„ 



V B 



(54) 



\b\/Va' 

where A is the cell are^ Our results are highly divergence- 
free, as we illustrate with steady-state solutions to three 
problems. For the Michel rotating monopole (Section 4.2.1 1, 



the normalised magnetic divergence is ~ 10 for an 84 x 



where Bang = {0, B^ B"^}, not including B' — we found that 56 grid, and ~ 10 for a gr id of 192 x 128. The twisted 



damping the radial magnetic field led to unphysical currents 
leaking back into the normal domain. This approach is sim- 
ilar to the matched-layer method of Yang et al. (1997| . The 
frictional coefficient, as, is chosen to be zero in most of the 
domain, and to rise smoothly near the boundary; the func- 
tional form we use is 







ctq 1 - exp 



^out 



if r < rs 
if r > 

(53) 



magnetosphere (Section 4.2.2 I has normalised V ■ B mostly 
around ~ 10~^^, rising to 10~^ in the region with largest 
current. Our fiducial aligned rotator solution (Section|5| has 
normalised V ■ B of ~ 10^* in and near current sheets and 
~ ]^Q-io elsewhere; the highest-resolution solution has values 
roughly ten times lower everywhere. 

The higher V • B near current sheets is at least partly 



Our discretisation is based on nodes rather than cells, but here 
it is unimportant which adjacent nodes are chosen to form a fic- 
titious cell. 
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due to the presence of the aforementioned Gibbs oscilla- 
tions, and would be expected to be much smaller in the 
recoverable high-accuracy solution (see Section |3.6[ ). This 
Gibbs-generated divergence would also be present if a vec- 
tor potential had been used. Finally, V • i? is higher than 
normal right at the inner boundary, presumably because of 
the filtering applied to the Chebyshev derivatives, which be- 
come increasingly dependent on the highest frequencies as 
the boundaries are approached (e.g. |Godon fc Shaviv|1993| ). 

Most importantly, no unusual or troubling behaviour 
has been observed to be correlated to an increase of the 
magnetic field divergence on the grid, the evolution appears 
to be stable for very long run times, and the measured di- 
vergence decreases with increasing resolution. 



3.10 Code infrastructure 

PHAEDRA is written in C for speed and portability. The 
spectral transforms are performed with the FFTW3 li- 
brary ( Frigo & Johnson 2005 1 , which allows transforms 



of arbitrary size with 0{N log N) complexity. The code 
is fully MPI-parallel, with a simple automatic domain- 
decomposition function which does not require the grid di- 
mensions to be a multiple of the number of processors. 

The parallelisation works similarly to the method of 
| |Pelz|199T| . The domain, in real space, is slab-decomposed 
in the radial direction. In the forward transform to spectral 
space, f{r,8) — >■ a„;, the ^-transforms are first performed 
using data local to each processor, after which the mixed 
fi{r) data undergo a parallel matrix transpose, done with 
a single MPI_Alltoallv() call. The r-transforms can then 
be performed, and the coefficients a^i, now decomposed in 
the Z-direction, used to calculate the coefficients of the de- 
sired derivative. The real-space derivative values are found 
by repeating the above steps in the reverse order. The filter 
is applied just before the inverse transform. In the SSV fil- 
tering step, performed at the end of each full time step, the 
coefficients are calculated, filtered, and immediately trans- 
formed back to real space. To date the code has been run 
on between one and 64 processors. 

The data output is done collectively, in parallel, us- 
ing the Parallel HDF5 library. The data are accessed via 
the XDMF standard, in which an auxiliary XML file, also 
written by the code, describes the contents of the HDF file 
containing the data arrays. This format allows the data to 
be easily opened by the Visit visualisation software, among 
others, without the need of a custom plugin. 



4 TEST PROBLEMS 
4.1 ID tests 

For the following two test cases we use a one-dimensional 
simplification of the code, with Chebyshev polynomials as 
the basis functions, an arcsine coordinate map, and a Carte- 
sian vector basis (which simply requires using the metric 



= Si 



The boundary conditions are enforced strongly 
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Figure 3. Stationary Alfven wave test problem. The line shows 
the exact solution, the dots show the numerical solution. 



4.1.1 Stationary Alfven wave 



Komissarov (2004a I describes an analytical solution for a 
stationary Alfven wave: Bx = By = Ez = 1, Ey — 0, 

(1 for X < 0, 

1 + 0.15 {1 + sin [57r {x - 0.1)]} for < x < 0.2, 
1.3 for X > 0.2, 

(55) 

and Ex ~ —Bz. For ease of comparison with this paper, we 
also use TV = 200 and a domain x G [-1.5,1.5]. The SSV 
filter strength is set to assv = 0.1 (even though no such 
filtering is required for this smooth problem) so that the 
effect of a discontinuity-capturing level of diffusion can be 
observed. 

Fig. [3] shows the numerical solution for Bz a,t t = 2 and 
t — 4000. The profile does not broaden noticeably, even after 
more than a million time steps. The small wiggles, on either 
side of the jump in the t — 4000 solution, are imprinted 
by the action of the SSV filter on a function which is not 
infinitely smooth {Bz has discontinuous first derivatives at 
j: = and X = 0.2). 



4.1.2 Riemann problem 

A Riemann problem which results in a current sheet is de- 
scribed by Komissarovj ( |2004a) . The initial conditions are: 
S = 0, - 



0, Bx = 1, and 

Byix) = 



Bo for a; < 0, 
Bo for a; > 0. 



(56) 



at n = 0, A'^ — 1; no sponge layers or non-refiecting boundary 
conditions are used. 



The current sheet forms spontaneously at a; = 0, and two 
fast step waves are emitted, one in either direction. The 
numerical solution for By at t = 1 is shown in Fig. [4] assv = 
0.1, and 2pssv = 8 as usual. The three jump discontinuities 
are clearly unresolved on the grid, and the fast waves remain 
sharp. The Gibbs oscillations are confined to the immediate 
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Figure 4. Riemann problem at t = 1, with current sheet at x = 0, 
for two grid sizes A^. 



vicinity of a discontinuity, and the affected region shrinks as 
resolution is increased. 



4.2 2D tests 

In these problems we used the full 2D code, in spherical 
coordinates. In all cases the star has a radius of unity, rin — 
1; its radial light-crossing time-scale is then At = 1. 



4-2.1 Michel monopole & split monopole 

The exact three-dimensional solution for the field configu- 
ration, in force-free electrodynamics, surrounding a rotating 



magnetic monopole was derived by Michel ( 1973 1. Given ini- 
tial conditions B'' = /o/r^ = B"^ = 0, E = 0, and sub- 
ject to uniform rotation with angular velocity Q, the steady- 
state solution has a toroidal magnetic field 



Bj, 



sin 9 



(57) 



This analytic solution is well suited to a multi- 
dimensional convergence test. We set the domain to be 
1 ^ r ^ 30 and vary the grid size N x L, holding A'' ~ 1.5L. 
No sponge layer is used, since monopole field lines imply 
radial outgoing waves, for which the characteristic bound- 
ary treatment is very effective (Section 3.8.21. The order of 
the aliasing-controlling filter is 2p = 36 for L ^ 32, and 
slightly lower for smaller grid sizes; no SSV filtering is ap- 
plied because the solution is smooth. The angular velocity 
is smoothly increased from zero to f2 = 0.1 in twenty radial 
light-crossing times (i.e. between t = and t = 20), and 
the solution is sampled at t = 100, on a surface of constant 
radius at r = 5. The fractional errors in B^, defined as 



fractional error — 



B 



analytic 





(58) 



Figure 5. Fractional errors in for a rotating monopole. Grid 
sizes are, from top, L = 16,20,24,32,40,56; 1.5L. N and 

L are the number of nodes in the radial and angular directions, 
respectively. 



are plotted in Fig. [5] The errors decrease approximately ex- 
ponentially with increasing resolution, reaching a level of 
roughly lO"^'^ with a grid size of 84 x 56. No effort was 
made to optimise the domain size, or the radius or time at 
which the numerical and analytic solutions are compared. 

A similar test can be performed for the split monopole, 
which simply involves reversing the sign of /o , and therefore 
of B^ and B'^ in the solution, across the equator. The dis- 
continuous magnetic field implies the existence of an equa- 
torial current sheet. This configuration allows us to test the 
behaviour of the code in the presence of a realistic disconti- 
nuity; we use SSV filtering with assv ~ 0.05. Fig. |6] shows 
the errors for the split monopole, for a problem otherwise 
identical to that described previously. 

In this case the convergence is not uniform, being faster 
farther from the discontinuity at the equator. Order unity 
errors remain immediately beside the current sheet, as de- 
scribed in Section |3.5[ but the solution converges strongly 
everywhere else, reaching a fractional error floor of 10~^ 
for the 224 x 155 grid. The large pointwise errors near the 
equator are not worrying, because they are due to the super- 
position of well-behaved, understood, and controlled Gibbs 
oscillations on top of an otherwise accurate solution, rather 
than uncontrolled numerical noise or error; see Section [3.6| 
for a discussion. 



4.2.2 Twisted dipole 

As another 2D test, consider a dipole magnetosphere that is 
twisted by a differential rotation of the star's surface. We as- 
sume a latitude-dependent surface motion that is symmetric 
about the dipole axis; then the magnetosphere remains ax- 
isymmetric. The twisted field lines become extended in the 
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Figure 6. Fractional errors in for a rotating split monopole. 
Grid sizes are, from top, L = 45, 65, 95, 155; N 1.5L. 
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r 

Figure 7. versus r at d = ■it/2, at t = 200, once the twist 
profile has been implanted; simulation (a) (see text). 



azimuthal direction between their footpoints in the northern 
and southern hemispheres. 

We start from a normal dipole configuration at t = 
and gradually impart a twist by moving the surface; at t = T 
the surface motion stops. If the magnetosphere is an ideal 
conductor, it must remain static sd, t > T — the twist will 
remain 'frozen' and persist forever. In the presence of re- 



sistivity, the magnetosphere must untwist with time (Be 
loborodov|2009 1. Our ideal force- free model has no physical 
dissipation, and the rate of untwisting at t > T will measure 
the numerical dissipation for our scheme. 

The speed of the surface motion is much smaller than 
the speed of hght, and so the twist will be implanted slowly, 
relative to the relevant wave-crossing time-scales. In our sim- 
ulation the stellar rotation is applied only in an annular re- 
gion of one hemisphere: 

"^^^ " i + e^p{-K[g{e- e,) + euw]}' ^^^^ 

where 9c is the colatitude of the centre of the annulus, 6hw 
is its angular half-width, k determines the sharpness of its 
edges, and 



9 ' 



1 ae <9c 
-1 if6l>6'c 



(60) 



The rotation rate is increased smoothly from zero to Sin 
in time T/2, and returned symmetrically to zero. 



(n„,ax/2) [1 - cos (27rt/r)] if t T, 
if t>T, 



(61) 



twisting the affected region of the star, — ^hw ^ ^ ^ 
Oc + Ohw, through an angle of -0 = ilniaxT/S. 

Here we use a grid of A'' x L = 320 x 255 on a do- 
main r £ [1,40], and twist an annular region of the north- 
ern hemisphere, given by 9c — O.Mvr, Onw = 0.047r, and 



K = 60. For this test we want only a small perturbation 
on top of a dipole field (small total twist amplitude), and 
for the magnetosphere to move slowly through a sequence 
of quasi-equilibrium states (requiring small flmax and large 
T), and so choose f2max ~ 10~^ and T — 200. This implants 
a twist of i/) = 0.1 radians by time t = 200, after which the 
stellar surface is at rest. 

To highlight the importance of the parallel current (see 



Sec. 3.7 1 we perform three simulations, all using the usual 
eighth-order SSV filtering with assv = 0.05. In run (a) we 
use the full force-free current in the equations of motion, J — 
J\\ + Jdrift from equation (40 1; in (b) we keep only the drift 
current, J = Jdrift, and remove accumulated at every 
Runge-Kutta substep; in (c) we use only the drift current, 
and remove parallel electric field at the end of every full time 
step. Run (a) therefore uses our standard numerical scheme; 
for this simulation the profile of along the equator at 
t = 200 is shown in Fig. [7| 

The total twist along any field line can be calculated 
by numerically integrating its path in space: the twist is 
tp ~ |0ond — <?!'start|. We label the field lines with the frac- 
tional fiux function, u. It measures the fraction of the star's 
magnetic flux passing through a circular contour, centred on 
the magnetic axis, which goes through the field line's foot- 
point. This function is m = sin'^ 6i for a dipole field, where 
01 is the co-latitude of the northern footpoint, a relation- 
ship which is unchanged by any axisymmetric motion of the 
stellar surface. 

Fig. [5] shows the measured twists at two times, t = 200 
(solid lines) and t = 4000 (dashed lines), for each of the 
three runs. In run (a) the curves lie on top of one another; 
in fact, sharp twist profiles are preserved, at this resolution, 
even for tens of thousands of light-crossing times. In run (b) 
the profile at t — 200 is fairly good, but it has diffused sig- 
nificantly hy t — 4000, while in run (c) the profile is already 
inaccurate by the earlier time and has almost completely dis- 



© 0000 RAS, MNRAS 000, 000-000 





logio * 



Figure 9. Energy injected into magnetosphere by twisting, in 
units of initial field energy: Wi — Wq (solid line), W2 (dashed 
line). The lines lie on top of one another. 



appeared by the later. These three cases demonstrate how 
important the parallel current is for self-consistently main- 
taining E ■ B — throughout each time step; omitting it 
implies making errors which demand the removal of elec- 



tric field from the system, introducing an artificial source of 
dissipation. 

We focus now on run (a). The permanence of the im- 
planted twist seen in Fig. |8] can also be interpreted in terms 
of energy conservation. Define Wi as the total electromag- 
netic energy of the twisted configuration, instantaneously 
measured on the grid, 

Wi^^J{B^ + E') AV, (62) 

where V is the three-dimensional volume of the computa- 
tional domain. Let Wo be the energy for the initial pure 
dipole field. Define W2 as the energy calculated by numeri- 
cally integrating the net energy flux into the grid over time, 

W2^ [ [Un{t') - Lout(t')] dt', (63) 
Jo 

where L = lo ^ ■^)'" sin^dS, and the subscripts de- 
note the luminosity measured at the inner and outer bound- 
ary surfaces respectively. The part of the grid containing the 
absorbing layer is not included in calculating Wi or 14^2- If 
there is no dissipation, we expect Wi — Wo = W2. 

The measured twist energy Wi — Wo and the im- 
parted energy W2 are shown in Fig. |9] as a function of 
time. The two curves lie on top of one another to within 
their widths. The energy on the grid is stable once the 
twist-up phase is complete. Between t — 200 and t = 
4000 the twist energy decreases by a factor of 4 x 10"*. 
At t — 4000, the fractional difference between accumu- 
lated twist energy and integrated net Poynting flux is 
[(Wi - Wo) - W2] I {Wx - Wo) = -3.9 X 10"*. 
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5 THE ALIGNED ROTATOR 



5.2 Evolution to steady state 



The objects which we observe as radio pulsars are gener- 
ally accepted to be rotating magnetised neutron stars, with 
magnetic fields of 10^^ G or more. If the star's basic mag- 
netic field is modelled as a dipole, the electric field induced 
by its rotation, equation ( |44[ ), has a large radial component, 
which is able to rip charge d particles from the stellar surface 
(Goldreich & Julian 1969|. These particles, and e* pairs cre- 
ated in the magnetosphere, surround the star with force-free 
plasma. 

To simplify the problem, one can consider the axisym- 
metric configuration with the magnetic moment parallel to 
the rotation axis: the aligned rotator. A steady-state so- 
lution was found by Contopoulos et al. (19991, which in- 



cluded a region of closed field lines extending to the light 
cylinder (defined as the cylindrical radius i?LC at which the 
co-rotation speed is the speed of light), and open, asymptot- 
ically monopolar, field lines extending to infinity. A current 
sheet is present at the equator beyond the light cylinder, 
which splits at the 'Y-point' to follow the last closed fiux 
surface. It was later found that equilibrium solutions exist 
with the Y-point at any distance from the rotation axis. 



within the light cylinder (Goodwin et al. 2004 Contopou- 



los 2005 Timokhin 20061. Absent a unique solution, one 



turns to time-dependent studies (Spitkovsky 


2006 


McKin- 


ney||2006b| |Komissarov|2006| |Kalapotharakos & Contopou- 


los 


2009 


Yu 


2011 


1, which have all found that the Y-point 



moves toward the light cylinder, where it subsequently re- 
mains. 



5.1 Numerical setup 

In our simulations, spatial and temporal scales are fixed by 
setting r* = 1, where r, is the stellar radius; c = 1 as 
usual. The star is smoothly spun up from rest, over a few 
light-crossing times, to the rotation frequency il, implying 
a light cylinder at i?LC = 1/^- We have investigated cases 
with _Rlc ~ 5, 10, 20, 30, and 40; there were only minor 
differences between the solutions, and here we concentrate 
on those with i?LC = 20. In these we set rout = 70, with 
an absorbing layer beginning at r = 50. A radial arcsine- 
plus-algebraic coordinate mapping is used, with parameters 
e — 10~^^ and Q = 0.7; the grid is equally spaced in the 
meridional direction. The current sheets are captured with 
super spectral viscosity (SSV): 2pssv = 8, assv = 0.1. 

We have performed simulations with a range of grid 
sizes, from A'' x L = 81 x 49 to 768 x 507. In order to demon- 
strate that a very fine grid is not needed for accurate results, 
we will concentrate on our run with a grid of 384 x 255 — all 
the following results are from this simulation unless explic- 
itly noted otherwise. The behaviour is similar for all resolu- 
tions. 

The initial magnetic field is set to a dipole with unit 
magnetic moment: Br — 2cos6/r^, Bo = sinS/r"^, B^ — 0. 
Rotation is introduced by applying an electric field, equa- 



During the spin-up phase, Alfven waves are injected along 
field lines into the magnetosphere, filling it with charges and 
currents. The magnetic energy increases and the poloidal 
field lines inflate, appearing to be pulled outward along the 
equator, where waves from opposite hemispheres meet. By 
t = 40 a clear distinction can be seen between those lines 
which are too close to the star to be strongly inflated, and 
those which are on the path to opening. Waves can be seen 
to propagate backwards and forwards on the former, while 
the latter are expanding at nearly the speed of light. Poloidal 
field line projections are drawn in Fig. |10[ 

Let us define the closed zone as that region with oscil- 
lating fields and currents, and the Y-point as the point at 
which the currents around the closed zone meet at the equa- 
tor. All the fields are smooth before t ~ 35, at which time 
the radial and azimuthal currents collapse to an unresolved 
equatorial current sheet. As the configuration evolves, all 
components of the current around the closed zone, and in 
the equatorial current sheet, change direction, as one can 
see in Fig. |11| The Y-point is at r = 0.5-Rlc before the 
reconfiguration, and at r = 0.6-Rlc afterwards. By t = 90 
Alfven waves have again filled the closed zone up to the Y- 
point, and the quasi-steady march of the Y-point to the light 
cylinder begins. 

The equatorial current sheet is unresolved on the grid, 
having discontinuities in both Br and B^; Fig.[2]^a) shows B^ 
against 9 at r = 22.25. The total magnetic field drops close 
to zero directly on the equator, and so the magnitudes of any 
electric fields must be reduced so that the second force-free 



tion (441, bringing the star to its final angular velocity by 
t = 10. 



condition, equation (131, is not violated (see Section 3.7 1. 
When choosing the grid, it is advantageous to have a line 
of nodes on the equator, otherwise the current sheet must 
choose which of the two equidistant nodes to collapse onto — 
since neither choice obeys the symmetry of the problem, the 
sheet will periodically move from one line to the other. 

At t ~ 100, more open flux has been created than can 
be supported by the energy being pumped into the mag- 
netosphere by the star's rotation. The Y-point moves very 
slowly towards the light cylinder, as some of the open fleld 
lines reconnect in the current sheet. The Alfven waves in 
the closed zo ne a re eventually damped by numerical dissi- 
patiorQ Fig. 12 shows the toroidal magnetic field for low 
(256 X 155) and high (768 x 507) resolution runs at t = 250. 
The oscillations survive longer with higher resolution, but 
the Y-point moves only slightly slower; here, it is at r = 16.8, 
versus 17.5 for the low-resolution run. 

The Y-point reaches the light cylinder at f ~ 500, but 
open fiux continues to reconnect slowly until t ~ 700, about 
5.5 rotation periods. After this time the solution is station- 
ary. The time taken to reach equilibrium is practically in- 
dependent of grid size. The equilibrium solution has some 
field lines closed outside the light cylinder due to the ef- 
fective resistivity of the current sheet; for example see the 
last panel of Fig. [lO] or Fig. [15] below. Note however that 
the zero-toroidal-field closed zone is strictly within the light 
cylinder. 



^ The waves are sheared by field line curvature, becoming longer 
and thinner. Numerical dissipation becomes significant when they 
approach the grid scale, and are attenuated by the filters. 
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Figure 10. Poloidal field lines for the aligned rotator, from non-rotating at t = to the equilibrium solution at t = 1000. The thicker, 
red field line closes at the Y-point in the equilibrium solution; the dashed line at i? = 20 indicates the light cylinder. There are 25 lines 
drawn from each pole, evenly spaced in colatitude between O.Olvr and O.IStt. 



5.3 Steady state 

We now concentrate on the obtained equilibrium solution. 
Unlike numerical solutions in previous works, ours does not 
exhibit plasmoid emission by the Y-point once it has ap- 
proached the Ught cylinder. The Y-point and current sheet 
are steady for long times; we ran a 256 x 155 simulation un- 
til t = 10, 000, or 80 rotational periods, without seeing any 
indications of Y-point instability. We found plasmoid emis- 
sion from the Y-point in only two circumstances. Firstly, if 
the level of filtering was too low; no plasmoids were seen 
if the filtering was strong enough to prevent the Gibbs os- 
cillations on either side of the current sheets from growing 
over time. Secondly, if the radial grid spacing was too large 



near the light cylinder; we found the Y-point was stable if 
Ar < 0.75 rA6. This is probably due to the action of the 
spectral filters, which make the current sheet mildly resis- 
tive, slightly diffusing the Y-point. 

The current sheet resistivity largely comes from filtering 
in the meridional direction, across the sheet. As the sheet 
forms, B^{6) goes to zero at the equator, with a cusp-like 
profile on each side, implying significant high-wavenumber 
content. The filters damp these high wavenumbers, causing 
the smoothed to be non-zero on the equator, which closes 
field lines. Eventually an equilibrium is reached between fil- 
tering and the electromagnetic forces trying to compress the 
current sheet. We verified this picture by using unfiltered 
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Figure 11. Formation of ttie Y-point. Thie colour shows radial current density; gold is positive, blue is negative. Charge flows into, and 
accumulates at, the Y-point, until the Y-point reaches the light cylinder. The light cylinder is at i? = 20. 



values of whenever the — > condition is vio- 
lated; this resulted in near-zero magnetic field in the cur- 
rent sheet, and a Y-point that moved outwards much more 
slowly. However, without filtering the evolution eventually 
became unstable. 

To investigate the dissipation in the current sheet, we 
performed a simulation using a 576 x 255 grid, with an 
outer boundary at r = 1000, an absorbing layer begin- 
ning at r = 700, and a more severe coordinate mapping; 
again, i?LC = 20. In Fig. [l3] we plot integrated Poynting 
flux through concentric spheres, as a function of radius. The 
outgoing flux is constant within the light cylinder, with max- 
imum fractional variation of about 3 x 10~* near the star. 
The flux inside the light cylinder agrees with the value found 
in previous works, fi^Q.'^/c^, where fj, is the star's magnetic 
moment, to fractional accuracy of 6 X 10"^ for our fiducial 
simulation and 2 x 10""^ for our high-resolution one. 

Outside the light cylinder, some of the outgoing fiux is 
lost in the current sheet due to its effective resistivity due 



to the filters; this deficit decreases with increased resolution. 
The energy loss is relatively large, because the resistivity is 
confined to the current sheet, which is kept sharp by the 
fully ideal surrounding magnetosphere. Solutions with global 
resistivity dissipate a smaller fraction of their luminosity in 
the current sheet ( Kalapotharakos et al.|20TT I . As shown in 
Fig. |14[ the open field lines are asymptotically radial. 

The magnetic field lines closing inside the Y-point form 
the closed zone, with no poloidal currents and zero toroidal 
magnetic field. The density of field lines is higher around the 
boundary of the closed zone. This can be seen in Fig. |15| 
where we plot field lines and contours of toroidal magnetic 
field. The closed field lines outside the light cylinder form 
cusps at the equatorial current sheet (we do not consider 
these lines to be part of the 'closed zone' because particles 
are not trapped on them, being able to flow out through the 
current sheet). The total magnetic flux through the light 
cylinder, in one hemisphere, is 1.39 ± 0.01 fiQ./c. This is 



larger than the value, 1.23, obtained by Contopoulos (20051 
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Figure 12. Comparing low resolution (256 X 155), left, and high resolution (768 X 507), right, simulations at t 
toroidal magnetic field. The Y-point has not yet reached the light cylinder at = 20. 



250. The colour is 




Figure 13. Luminosity (integrated Poynting flux) through a 
sphere of radius r, as a function of r. is the luminosity mea- 
sured at the stellar surface; it equals fj^W^/c^ (to a fractional 
accuracy of 6 x 10~^). 



and Timokhin ( 2006 1 for ideal steady-state force- free mod- 
els, because some of the closed fiux in our solution has dif- 
fused through the light cylinder due to finite resistivity. 

The magnetosphere has current leaving the polar cap 
and, outside the Hght cylinder, returning to the star mostly 
in the current sheet. At the Y-point the current sheet splits 
in two, and follows the boundary of the closed zone. In 
Fig. [it] we plot the poloidal current density measured on 




200 300 



Figure 14. Poloidal fleld lines, out to r = 35i?LC- The field lines 
have the same footpoints as in Fig. |10| 



the star, normalised to the Goldreich- Julian current density 
(the speed of hght times the equilibrium charge density); a 
similar plot was obtained for the steady-state solution by 



Timokhin (20061. To aid comparison with Fig. 5 of that pa- 



per, we have scaled the axis to the width of the polar cap. 
The current sheet is seen just inside 9/6pc = 1. 

The three components of the current density, and the 
charge density, of the equilibrium solution are shown to a 
common scale in Fig. |16| Some Gibbs oscillations are ap- 
parent near the Y-point and beside the outer current sheet, 
but they are controlled by the filtering and do not appear 
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Figure 15. Poloidal field lines and contours of constant B^. 
There are 100 lines drawn from each pole, equally spaced in co- 
latitude between 6 = O.OOStt and O.IItt. Contours, 14 in each 
hemisphere, are drawn in blue for negative values and red for 
positive values. They show IB^I equally spaced between 3 X 10~^ 
and 10-3. 



to cause any problems. The closed zone has no poloidal cur- 
rent, but there is toroidal current from the corotating charge 
density. 

In a steady, ideal, force-free magnetosphere the elec- 
tric field has a potential that is constant along the magnetic 
field lines. Our solution is everywhere close to this behaviour 
(except in the equatorial current sheet where numerical re- 
sistivity is significant). The equipotentials follow the shape 
of magnetic field lines, and are smooth and accurate. The 
charge density corresponds to the Laplacian of the poten- 
tial, and the second derivatives have more numerical noise. 
Nevertheless, the obtained charge density reproduces all the 
expected features, including steep gradients near the current 
sheet and the singular behaviour at the Y-point. The cur- 
rent sheet is positively charged outside the Y-point and neg- 
atively charged around the closed zone, in agreement with 



the previously obtained steady-state solutions (see Timo- 
|khin| [2006| . The negatively-charged current sheet around 
the closed zone appears to be resolved, and the thickness 
of the negatively-charged region decreases slowly with in- 
creased resolution. This suggests a thickening of the current 



sheet due to finite resistivity, as argued by Gruzinov (2011 1 



The thickening must occur due to resistivity near and out- 
side the Y-point, as dissipation inside the hght cylinder is 



significant voltage. In most observed pulsars, this voltage is 
not so large as to make the force-free approximation unrea- 
sonable. A real danger for the force-free model appears if 
the required charge density or current cannot be created. 
The electron-positron dischar ge operates at the polar 



cap if a ^ < 1, where a = J/cpe ( Beloborodov|[2008 1 . This 



condition is satisfied in the zone of return current (J^ > 



0), near the edge of the polar cap where a < (Fig. 171. 
Pairs are created with a high multiplicity and outflow along 
the magnetic fleld lines, screening E ■ B. The presence of 
dense pair plasma makes the force-free approximation safe 
in the return-current zone (except in the current sheet). The 
boundary of this zone (Jr = 0) is shown by the blue curve 
in Fig. [1|] 

The discharge does not occur in the central parts of the 
polar cap where Jr < and < a < 1. Instead, the required 
pe and J are supplied by the charge-separated flow pulled 
out from the star (Beloborodov 2008; Chen & Beloborodov, 
in preparation) . The force- free approximation remains accu- 
rate along the field lines extending from this region as long 
as < Of < 1 remains satisfied. We find that a < 1 every- 
where in the zone Jr < 0, inside and outside of the light 
cylinder. However, a > is not satisfied. There is a small 
region in the zone Jr > outside the hght cylinder where 
a becomes negative (because pe changes sign, see Fig. 18 1. 



The charge-separated outflow passing through this region 
fails to supply the charge density of the required sign, and 
a large E ■ B must develop. E ■ B may be limited by locally 
initiated pair creation in young, fast-rotating pulsars; how- 
ever, in most pulsars pair creation is inefficient so far from 
the neutron star. 

We conclude that, for the aligned rotator, the force- free 
approximation is expected to fail in the region where Jr < 
and pe > 0. This region is, however, small, and this prob- 
lem may not impact the obtained global solution. We note 
that a similar, but larger, region is seen in Fig. 4 of lrTmokhinl 
(2006[) and in Fig. 3 of|Contopoulos et al.|([l999l). The differ- 



ence between their and our solutions is due to the fact that 
their model is strictly ideal everywhere, while our model is 
(nearly) ideal only outside the equatorial current sheet. The 
dissipation in the current sheet affects the magnetosphere 
as discussed above and shrinks the region that is dangerous 
for the force-free model outside the current sheet. 

We note also that the polar-cap accelerator in the zone 
of return current Jr > can supply some e^ pairs to the 
neighbouring field lines with Jr < (near the blue curve 
in Fig. 18 1. Pair creation is not exactly local to the accel- 



negligible (see Fig. 13 1 



eration region because it involves an intermediate step-the 
emission of a high-energy photon which must propagate a 
finite distance across the field lines before converting to e* . 
Pairs created on, and outflowing along, the field lines slightly 
outside the return-current zone may supply enough positive 
charges to the problematic small region (Jr < 0, > 0) 
and validate the force- free condition there. 



5.4 Viability of the force-free model 

The force-free model relies on the availability of charges that 
sustain the required charge density and electric currents. 
Charges can be pulled out from the star or supplied by e* 
pair creation. Both mechanisms require a longitudinal volt- 
age, i.e. E ■ B 0; pair creation, in particular, requires a 



6 DISCUSSION 
6.1 PHAEDRA 

We have found that phaedra converges exponentially 
quickly to smooth solutions, is stable and accurate when 
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discontinuities are present, and exhibits very low numerical 
diffusion outside of current sheets. A critical ingredient is 
the full force-free current density, in particular the parallel 
conduction current, which maintains E ■ B = without in- 
troducing any stiffness to the equations of motion. We are 
able to use the full current because our mesh is unstaggered, 
and so its evaluation doesn't require any interpolation of 
fields or their derivatives. The unstaggered grid is itself en- 
abled by the inherent accuracy of the pseudospectral spatial 
derivatives. 

The flexibility of the mapped-coordinate method allows 



efficient calculation in large domains, while retaining accu- 
racy near the stellar surface at the inner boundary. It also 
permits resolution to be concentrated where large gradients 
are expected to form, a kind of fixed mesh refinement. With 
the adjustable parameter in the arcsine mapping of the ra- 
dial Chebyshev series, a deliberate tradeoff can be made 
between accuracy close to the boundary and stable time 
step. The ability to model interactions of force-free waves 
with a solid boundary to high accuracy will be useful when 
studying, for example, turbulence driven by neutron star vi- 
brations. The treatment of non-periodic boundaries is one 
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Figure 17. Normalised poloidal current density across the polar 
cap. Jq] = pGJi PGJ = — 2f2 ■ B, where c = 1 and 4it = 1. 
^pc = 0. 088-71 is the half-width of the polar cap, defined by the 
footpoint of the last field line in the closed zone. 




Figure 18. The blue line separates the outflowing current region, 
with Jr < 0, near the poles, and the return current region nearer 
the equator. The orange line is the contour of zero charge density 
pe, and the shaded region has pe > 0. Poloidal field lines are 
drawn in black. The charge density has been cleaned with the 
DTV filter (Section |3.6| to remove Gibbs oscillations near the 
equatorial current sheet. 



of the advantages of spectral expansions over high-order fi- 
nite difTerence methods, which can experience difficulties at 
surfaces. 

The code is parallelised with the standard MPI library, 
and can be run in both shared- and distributed-memory en- 
vironments. It is efficient enough to run on dense grids, using 
around 16-32 processors, for many millions of time steps. 
Simulations on coarser grids can be run on one or a few 
processors, on consumer laptops and workstations. One con- 
cern is that the 0{N^) communication time required for the 
global MPI call will become a problem when scaling to a 
very large number of processors, for example for a three- 
dimensional calculation. We expect that the code should 
scale well up to several hundred processors on existing hard- 
ware, and that this will be sufficient. 

The principal issue we have encountered is the resistiv- 
ity imparted to current sheets by the spectral filters. There 
are numerous tricks that can be used to reduce or elimi- 
nate this in simple cases, like the aligned rotator we de- 
scribe above. However, we prefer not to use any method 
insufficiently robust, flexible, or efficient to also be applica- 
ble to general three-dimensional current sheets propagating 
across the grid. The ideal aligned rotator may represent the 
worst-case comparison of PHAEDRA to codes employing 
finite-difference or finite-volume methods, which can evolve 
the solution to a steady-state with very little flux closed out- 
side the light cylinder, either with use of a staggered mesh 



I Spitkovsky 2006 1 , or manual nulling of the inflow of flux 



into the current sheet (McKinney 2006b I. We are investl 



gating adaptive spectral- and real-space filtering, which may 
be able to stabilise current sheets with less dissipation than 
global filtering. 

In any case, physical current sheets are believed to pos- 
sess finite resistivity (e.g. Zenitani fc Hoshino|2007 Zenitani 
fc Hesse]|2008 1. The effective hyper-resistivity of the filters 
confines the dissipation to grid-scale features like current 
sheets where it is physically expected, leaving resolved fea- 
tures to evolve ideally. In this sense, our aligned rotator 
solution may be realistic. 



6.2 Planned Extensions 

We are planning to extend the scheme in several ways. 
The derivatives are already being calculated in a metric- 
independent fashion, and so it should be straightforward to 
adapt the code to the Schwarzschild or Kerr metrics, using 



the formalism of Komissarov ( 2004aJ . Eventually it should 
even be possible to include effects due to a time-dependent 
metric in this framework ( Komissarov|201~ I, possibly using 
data from a code which evolves the Einstein equations. 

The method currently assumes axisymmetry; we intend 
to make it fully three-dimensional, using either the complete 
'double Fourier' method on spheres or expansions in spher- 
ical harmonics. We expect the overall scheme to be adapt- 
able to different choices of basis functions with only minor 
changes. This would allow us to investigate the oblique pul- 
sar magnetosphere in spherical coordinates, and so alleviate 
some of the difficulties, such as stair-stepping on the in- 
ner surface and an inflexible equispaced Cartesian grid, en- 
countered by some existing 3D force-free codes. Removing 
the axisymmetric restriction will also permit us to evalu- 
ate the stability of field configurations produced by general 
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surface footpoint motions. Considering the lower resolution 
that is possible with 3D calculations, it is encouraging that 
we see close agreement between simulations on coarse and 
fine grids, even for the pulsar solution. 

Another promising modification is to a 3D Cartesian ge- 
ometry, using Fourier series in all directions. This geometry 
would be suitable for studying ultra-relativistic turbulent 
cascades ( Thompson fc Blaes|1998 Cho|2005 1 and instabili- 
ties surrounding nearly force-free current sheets. Aside from 
being simpler than expansions in spherical coordinates, the 
Cartesian-plus-Fourier combination has the benefit of allow- 
ing the solenoidal constraint on the magnetic field to be eas- 
ily enforced in spectral space. 

The low intrinsic numerical diffusivity of our code makes 
it ideal for studying the effect of physical plasma resistiv- 
ity. Although force-free electrodynamics lacks a well-defined 
fluid frame, and hence a preferred form for any dissipa- 
tive terms, several formulations of resistive, nearly force- free. 



electrodynamics have been proposed ( Lyutikov|2003| Gruzi- 



nov|2008[[Li et al.|201l"] [Kalapotharakos et al.|201li . Being 
able to resolve waves with many fewer points per wave- 
length, spectral codes require less diffusivity than lower- 
order schemes to accurately transport, without unphysi- 
cal oscillations, a given profile on the same grid. |Brandeii^ 
burg ( 2003 1 found that spectral derivatives permitted a vis- 



cosity five times lower than that needed for even a sixth- 
order finite-difference method. If enough physical resistivity 
is added to resolve otherwise sharp current layers, we will 
probably be able to dispense with the 'super spectral viscos- 
ity' filtering, and the code should again achieve exponential 
convergence. 
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